/*
  This function calculates 1st derivative, du/dr, with "ndiscr"-th order accuracy. It accepts "nxpt" number of
  grid-points co-ordinates "x" and "nxdervpt" number of derivative-points "xdervpt", where the 1st derivate
  is to be calculated. It returns the co-efficient matrix "mtrx" of size [nxdervpt]x[nxpt]. It doesn't filter
  out the zero-elements. */

# include <stdlib.h>
# include <math.h>

static void lagrange_derv(double *lderv, double *xstn, unsigned short nstn, short ipstn, double dxstn);

/* This evaluates the co-efficient matrix for ndiscr-th oder discratization of the 1st-derivative. */
extern short cmtrx_CD(double *x, unsigned short nxpt, double *xderv, unsigned short nxdervpt,\
 unsigned short ndiscr, double **mtrx)
  {
   unsigned short j, *tmtrx;
   short k, idxmin;
   double dx, dxmin, epsilon;

/* The discretization order has to be even. */
   if(ndiscr%2 != 0) return -10;

/* This the transformation matrix and it is used to mark the beginning of the stencil
   on which the discratization will be carried out. It becomes very useful when we
   approach towards ends of the domain. */
   tmtrx = (unsigned short *)malloc(sizeof(unsigned short)*nxdervpt);

/* This is used to check if two float quantities are same or not. */
   for(epsilon = fabs(x[nxpt-1]-x[0]), j = 1; j < nxpt; j++) if(epsilon > (dx = fabs(x[j]-x[j-1])) ) epsilon = dx;
   epsilon /= 100;

   for(k = 0; k < nxdervpt; k++)
      {
    /* Finding out the closest grid-point from the k-th derivative point. */
       dxmin = x[nxpt-1]-x[0];
       idxmin = -1;
       for(j = 0; j < nxpt; j++)
          {
           mtrx[k][j] = (double)0; /* Initializing the matrix's k-th row with zero. */

           if(fabs(dx = xderv[k]-x[j]) < fabs(dxmin))
              {/* Locating the closest grid-point from the k-th derivative-point. */
               dxmin = dx; /* The closest distance. */
               idxmin = j; /* The grid-number corresponding to the closest grid. */
              }
           }

   /* Check if the k-th derivative-point coincides with any of the grid-points. */
      if(fabs(dxmin) > epsilon)
          {/* The k-th derivative-point doesn't coincide with a grid-point. */
           if(nxpt < ndiscr) return -11;

        /* Get the grid-number right before the k-th derivative point. */
           if(dxmin < 0) idxmin--;

        /* Mark the beginning of the stencil that will be used in evaluation of the derivative. */
        /* The k-th derivative point is close to beginning of the grids. */
           if(idxmin <= ndiscr/2-1) tmtrx[k] = 0;
        /* The k-th derivative point is close to end of the grids. */
           else if(idxmin >= nxpt-ndiscr/2-1) tmtrx[k] = nxpt-ndiscr;
        /* The k-th derivative point lies in middle of the grids. */
           else tmtrx[k] = idxmin - ndiscr/2 + 1;

        /* This evaluates the derivative using the Lagrange polynomial at the k-th derivative-point. */
           lagrange_derv(&mtrx[k][tmtrx[k]], &x[tmtrx[k]], ndiscr, -1, xderv[k]-x[tmtrx[k]]);
          }
      else{/* The k-th derivative-point coincides with a grid-point. */
           if(nxpt < ndiscr+1) return -12;

        /* Mark the beginning of the stencil that will be used in evaluation of the derivative. */
        /* The k-th derivative point is close to beginning of the grids. */
           if(idxmin <= ndiscr/2) tmtrx[k] = 0;
        /* The k-th derivative point is close to end of the grids. */
           else if(idxmin >= nxpt-ndiscr/2-1) tmtrx[k] = nxpt-ndiscr-1;
        /* The k-th derivative point lies in middle of the grids. */
           else tmtrx[k] = idxmin-ndiscr/2;

        /* This evaluates the derivative using the Lagrange polynomial at the k-th derivative-point. */
           lagrange_derv(&mtrx[k][tmtrx[k]], &x[tmtrx[k]], ndiscr+1, k-tmtrx[k], 0);
          }
     }

   free(tmtrx);

   return 10;
  }



/* This function calculates a row-vector of the co-efficient matrix responsible for evaluation of the
  1st order derivative at a given point. It usage the Lagrange-polynomial approach and the row-vector
  consists of derivative of the Lagrange-weights at each stancil-point.

  xstn: Beginning of the stencil-points.
  nstn: Number of stencil-points. i.e. stencil-size.
  ipstn: If the point at which the derivative is evaluated coincides with grid-points,
         it marks that grid-point on the stencil. "dxstn" is ignored in this case.
  dxstn: If the derivative-point doesn't coincide with grid-points, it marks the distance
         of the point on the stencil. ipstn = -1.

  "ipstn" and "dxstn", both are relative to beginning of the stancil. */

static void lagrange_derv(double *lderv, double *xstn, unsigned short nstn, short ipstn, double dxstn)
  {
   unsigned short i, j;
   double wt, xp;

/* Locating the point on the grid at which the derivative is to be calculated. */
   xp = (ipstn < 0 ? xstn[0]+dxstn : xstn[ipstn]);

/* Itterating through all stencil points. */
   for(j = 0; j < nstn; j++)
      {
       wt = (double)1;

    /* Calculating weight of the Lagrange polynomial if the derivative-point doesn't conside with a grid-point. */
       if(ipstn < 0) { for(i = 0; i < nstn; i++) if(i != j) wt *= (xp-xstn[i])/(xstn[j]-xstn[i]); }

    /* If the derivative-point is a stencil-point but doesn't lie on this particular j-th stencil-point. */
       else if(ipstn != j)
          {
           lderv[j] = (double)1/(xstn[j]-xp);
           for(i = 0; i < nstn; i++) if((i != ipstn) && (i != j)) lderv[j] *= (xp-xstn[i])/(xstn[j]-xstn[i]);
           continue;
          }

       for(lderv[j] = 0, i = 0; i < nstn; i++) if(i != j) lderv[j] += wt/(xp-xstn[i]);
      }
   return;
  }
